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We investigate the disorder dependence of the static density, amplitude and current correlations 
within the attractive Hubbard model supplemented with on-site disorder. It is found that strong 
disorder favors a decoupling of density and amplitude correlations due to the formation of supercon¬ 
ducting islands. This emergent granularity also induces an enhancement of the density correlations 
on the SC islands whereas amplitude fluctuations are most pronounced in the ’insulating’ regions. 
While density and amplitude correlations are short-ranged at strong disorder we show that current 
correlations have a long-range tail due to the formation of percolative current paths in agreement 
with the constant behavior expected from the analysis of one-dimensional models. 
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I. INTRODUCTION 

More than 50 years ago Anderson has discussed the 
behavior of a superconductor in the presence of strong 
disorder.li] According to his analysis (and under the re¬ 
striction to elastic scattering from non-magnetic impuri¬ 
ties) the BCS wave-function, build from Bloch-type wave 
functions with opposite momenta, can be generalized to 
pairs made from the exact single-particle wave-functions 
of the disordered system plus their time-reversed partner. 
As a result one would expect a gradual dependence of the 
superconducting transition temperature on the presence 
of non-magnetic impurities caused mainly by a modifi¬ 
cation of single-electron properties as density of states 
etc. While this picture is certainly correct for weak dis¬ 
order, experiments on thin films of strongly disordered 
superconductorPlti^ have revealed a much more interest¬ 
ing behavior than suggested by Ref. [1] In particular, 
the observation of a superconductor-insulator transition 
(SIT) with increasing disorder provides evidence for an 
interesting interplay between localization of Cooper pairs 
and long-range superconducting (SC) order.E^ Moreover, 
the observation of a pseudogap in strongly disordered SC 
fjlm^UEl bears some resemblance to similar experimental 
findings in high-temperature superconductors'^"!^ that 
may suggest a common mechanism in some regions of 
the phase diagram. 

Theoretical investigations of disordered superconduc¬ 
tors are either based on bosonic or fermionic ap¬ 
proaches. In case of s-wave superconductivity the 
latter typically start from attractive Hubbard models 
where disorder is usually implemented via a shift of on¬ 
site energy levels These hamiltonians then are ei¬ 

ther treated within a sta ndard Bogoljubov-de Gennes 

approximatiorfi^ISlHlSIllIlel mo re soph isticated 

approaches like Monte-Carlo methods. h^B0 | 24 | gosonic 
models are then obtained from a large-U expansion, as 
e.g. the pseudospin XY model in a transverse fielcP^, 
where the hopping of Cooper pairs (corresponding to 
pseudospins aligned in the XY plane) competes with lo¬ 


calization due to random fields (corresponding to pseu¬ 
dospins aligned in the perpendicular direction). Further 
simplifications, as e.g. an Ising model in a random trans¬ 
verse field, are also introduced since they allow for ana¬ 
lytical treatments.!^ 

In recent years both approaches have lead to a coherent 
picture of the SIT: With increasing disorder the system 
starts to break up into “puddles” with finite SC order pa¬ 
rameter |A| > 0 and intermediate regions with |A| « 0 
although the spectral gap remains finite. The order pa¬ 
rameter distribution shows a universal scaling behavior, 
in agreement with experiment, where the relevant scal¬ 
ing variable is the logarithm of the order parameter dis¬ 
tribution normalized to its variance! ^^ * ^^ ! The phases of 
different puddles are weakly coupled, so that the system 
bears some resemblance with a granular superconductor. 
Upon applying a vector potential the system accommo¬ 
dates the phase twist in the regions with |A| « 0 so 
that the associated energy, and thus the superfluid stiff¬ 
ness, are strongly reduced. Moreover, calculations within 
the BdC approach of the attractive Hubbard model have 
shown that the induced current flows along a quasi one¬ 
dimensional percolative path or “superconducting back¬ 
bone” which connects the puddles.!^ This result has its 
counterpart in the analysis of the bosonic approach which 
has revealed a regime of broken-replica symmetry where 
the partition function is determined by a small number 
of paths.l^For both, fermionic and bosonic models, there 
exists a critical value for the disorder strength above 
which the system becomes insulating. The SIT is charac¬ 
terized by a vanishing of the superfluid stiffness, however, 
the single-particle gap persists across the transitiorP^. 

A still open issue is the nature of the spatial cor¬ 
relations in such granular SC state arising near the 
SIT. In the classical Ginzburg-Landau-Abrikosov-Gorkov 
theorjP^ there is a single scale ~ ^f/^, whose re¬ 
duction by disorder is mainly governed by the mean-free 
path i via ^ ~ On the other hand, in the vicinity of 

the Anderson localization transition the coher ence l ength 
is also controlled by the localization length. Con- 
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cerning the disordered attractive Hubbard model with 
a fragmented SC ground state as mentioned above, there 
is only limited knowledge about amplitude, density and 
current correlations. Previous Quantum Monte-Carlo 
studie^^ yield only limited information on the spatial 
dependence of the correlations due to the small (8 x 8) 
lattice sizes. On the other hand investigations of re¬ 
sponse functions on larger clusters within the BdG ap¬ 
proach where so far restricted to mean-field studies. In 
the present paper we evaluate the density, amplitude 
and current correlations by including fluctuations on top 
of the BdG solution thus generalizing the approach of 
Refs. I86I87I to the case with disorder. In particular we 
are interested in the question of how the physics is gov¬ 
erned by different length scales in different channels and 
how the formation of SC islands for strong disorder re¬ 
flects in the corresponding correlation lengths. 

The paper is organized as follows: The model is intro¬ 
duced in Sec. |ll] where we also outline the computation 
of correlation functions on the basis of the BdG ground 
state. Results are presented in Sec. [ITT] for amplitude, 
density and current correlations. We hnally conclude our 
discussion in Sec. m 




, compute the new values (T = 0) 


Ai = \U\y^Ui{n)v* (n) 

n 

K) = 2^|u,(n)p 


( 4 ) 

( 5 ) 


and iterate the obtained values, say K, (including also 
the chemical potential) up to a given accuracy 5K/K < e, 
typically e = 10“®. For the disordered systems studied 
in Sec. m clusters with up to 24 x 24 sites have been 
diagonalized. We mostly show results with filling n = 
0.875, but in some cases we also discuss the outcomes 
for smaller filling in order to avoid the proximity to half¬ 
filling, where specific effects can arise due to the tendency 
of the system to form a charge-density-wave (CDW) state 
as well. 


B. Amplitude and Charge Correlations 

We denote correlation functions by 

^ Idte^-\rdnm^m (6) 


II. FORMALISM 
A. BdG equations 

Our starting point is the attractive Hubbard model 
with local disorder 

H = ^ ) ^ijcla-'^jcr ~ \U\ ^ ^ ^ ( ViTlicr ( 1 ) 

ija i ia 

which we solve in mean-held using the BdG transforma¬ 
tion 

k 


where in the following 0,R correspond to either ampli¬ 
tude SAi or density Spi huctuations 

6Ai = {St]^ + St]1)/V2 

Spi = ^ ^ ^c|g.Cio. — (cJg.Cicr)^ , 
cr 

and we have dehned the pair huctuation operators 

^4 = 4^4^ - 

It is then convenient to dehne 2x2 matrices for the 
bare mean-held susceptibility 


UJkUn{k) = '^tnjUj{k) + [Vn - («n) “ p]Un{k) 

3 

+ A„u„(fc) (2) 

i^kVn{k) = - ^t^jVjik) - [Vn - - /r]u„(A:) 

3 

+ A>„(A:). (3) 

For simplicity only nearest-neighbor hopping tij = —t 
is considered in this work. The disorder variables Vi are 
taken from a hat, normalized distribution ranging from 
— Vq to -t-Vg. 

In the following Ui{k) and Vi{k) are taken to be real. 
Starting from an initial distribution of the gap Aj and 
density (rii) values we diagonalize the system of equations 




Azj Atj 
pA pp 
ij Aij 


( 7 ) 


and the interaction 


V- -\u\ 0 

' 0 -\U\/2 


( 8 ) 


which can be combined into “large” matrices according 


to 


/ Xll Xl 2 ’'' XlN \ 
0 0 D 

X2I X22 ■■■ X2N 

V X^l Xn 2 ' ■ ■ X%N / 


( 9 ) 
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and 


/V 0 ••• 0 \ 

0 y ••• 0 


V Q Q ■■■ V / 


( 10 ) 


The RPA resummation can then be written as 


x = x° + x°Zx 


which is solved by 


X = 


i-x°Z 


X 


( 11 ) 


Note that in this paper we will focus on static correla¬ 
tions. Since at Gaussian level the coupling between the 
phase fluctuations and the density/amplitude ones is pro¬ 
portional to the frequencjEH, they are decoupled in the 
static limit. On the other hand, the phase fluctuations 
enter in a crucial way in the calculation of the current 
fluctuations, as will be outlined in the next subsection. 


C. Current Correlations 


The current response J“(a;) to a vector potential 
Ax (n, uj) (which we fix along the a;-direction of our square 
lattice) is the sum of the diamagnetic and paramagnetic 
contributioiP^ 

H + x(i“0™)] ^x(to) (12) 


where G(n) = -t’E^icl ^.Cn+xcr + 4< 0 de¬ 
notes the kinetic energy on the bond between sites i?„ 

and i?„ -I- ax and j“ = Z)<T(cl+a.crCn<T - h.c.) is 
the operator of the paramagnetic current flowing from 
site Bn to Rn+a- Note that the notation for the cur¬ 
rent correlation function xUn^Jm) is slightly different 
from the correlations defined in the previous subsec¬ 
tion. At frequency w = 0 the current only couples to 
phase fluctuations = i{ST]i — Sr]l)/V2 via the vertices 
Km = X°(j“,<5$„) and Anm = Thus, the 

full (gauge invariant) current correlation function is then 
obtained from 


xUuJi) = x°(/“,Jm) 


V 


mk 


l-xK 


-1 


A 


•J kl 


Im 5 


(13) 


with x^ in the second term denoting the bare phase-phase 
correlation function and V/nfc = —USmk- 

For the Fourier transform of the configurational aver¬ 
age one finally obtains 


J“ = -D^^^Axiq) (14) 

where = -{Tx)Sa,x-{XqU°‘,r))- For J“ = and 
taking q along the y direction the limit limg^_>o = 


Dg corresponds to the superfluid stiffness and coincides 
with the quantity evaluated in Ref. [25] from an expansion 
of the mean-field free energy up to quadratic order in the 
vector potential. 



FIG. 1: (Color online) Top to bottom: amplitude (x^^(q)): 
density (x^^(‘l))i and off-diagonal (x'^^(q)) correlation func¬ 
tions in the superconducting state for parameter \U\/t = 2 
and in the clean limit (Fo/t = 0). 


III. RESULTS 

A. Correlations in the homogeneous system 

We start our considerations by a brief resume of the 
homogeneous case for which amplitude and density cor¬ 
relations have been analyzed in Ref. [33 and which are 
in agreement with our following finite cluster analysis. 
Fig. ^ shows the amplitude x^'^(*l)j density x^^(q) and 
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^pp 

(Q) at n = 0.875 is a remainder of this CDW insta¬ 
bility at half-filling which is transfered to the amplitude 
correlations via the mixed susceptibility shown in 

the bottom panel of Fig.Increasing \U\/t enhances the 
CDW correlations so that at some point the q = Q ampli¬ 
tude correlations also dominate with respect to the q = 0 
response. On the other hand the CDW correlations are 
suppressed in the dilute limit (not shown) so that upon 
reducing filling the maximum density response is first 
shifted away from Q along the Brillouin zone boundary 
and finally, below some concentration and depending on 
the value of |D|/t, the q = 0 response starts to dominate 
. A more detailed discussion on the filling dependence of 
the amplitude and density response in the clean case can 
be found in Ref. ITfl 


B. Disordered system: Real space analysis 


FIG. 2: (Color online) Distribution of the superconducting 
gap parameter Ai (displayed on a linear scale by circles) and 
superconducting currents (arrows) computed from Eq. (121 
for constant vector potential and a specific disorder con¬ 
figuration. Parameters \U\/t = 5, Vo/t = 2. The dashed line 
(blue arrow) indicates the cut which is analyzed in Fig. 


mixed correlation function for filling n = 0.875 

and \U\/t = 2 without disorder. 

For these parameters the maximum of the amplitude 
correlations is at q = 0 where it can be approximated as 


X^^(q) 


1 

m? -b cq^ 


(15) 


with the mass m and a parameter c characterizing the 
dispersion of excitations. The quantity = \/ c/w? can 
then be interpreted as a length scale for the decay of the 
amplitude correlations. On the other hand the density 
response is dominated by the contribution at q = Q = 
(tt, tt) and around this wave-vector can be described by 


x""(q«Q) 


1 

"t-q + CQ(q- Q )2 ' 


(16) 


In real space this corresponds to a staggered decay of the 
density correlations with length scale ^ cq/ttiq. 

The mixed susceptibility x^’^(q) is negative (positive) 
for densities n < 1 (n > 1) since the anomalous cor¬ 
relations {ci^Ci-^) are negative with a maximum of their 
absolute value at half-filling. Therefore a positive fluctu¬ 
ation in density 5p for n < 1 will lower (i.e. enhance the 
magnitude) the anomalous correlations. 

For the present model, in the absence of disorder 
and at half-filling, there is an “accidental” symmetrjESl 
which allows the superconducting order to be continu¬ 
ously rotated into the charge density wave (CDW) order 
at q = Q without energy change, promoting the charge 
density mode to a Goldstone mode. The enhancement of 


1. Mean-field solution 

For sizeable disorder the density varies on the scale of 
the lattice constant and correlates with the strongly spa¬ 
tially fluctuating disorder potential. Further on, it has 
been shown in Refs. ]m2i\22\2f,\2<^ that for strong disor¬ 
der the system disaggregates into SC islands with sizeable 
SC gap Ai which are embedded in regions with « 0. 
Fig.j^shows a map of the order parameter encoded on the 
size of the red circles showing the formation of the super¬ 
conducting islands. This island structure leads to a very 
weak superfluid stiffness. Indeed, upon applying a trans¬ 
verse vector potential, as done in Ref. [551 the current 
flows through an optimum percolative path or “super¬ 
conducting backbone” which determines the global stiff¬ 
ness. The latter not only depends on the volume fraction 
of the superconducting island, but also on the connec¬ 
tivity of these islands to the superconducting backbone. 
Thus one may have a moderate superconducting fraction 
and a very small global stiffness if the connectivity is 
poor. Fig. shows an example of the superconducting 
backbone for current circulation. Notice that it does not 
necessarily involve all significantly superconducting sites. 
For example, sites (1,7) and (3,12) in Fig. where 
is large, are left out which therefore are examples for 
poorly connected islands. Whereas connected islands de¬ 
termine the superfluid stiffness the disconnected islands 
dominantly contribute to the subgap absortpion in the 
optical conductivity.!^ 

Analyzing the mean-field solutions for several config¬ 
urations of disorder we find that there is a strong ten¬ 
dency to form superconducting dimers. For example, 
tor Vo/t = 2 ^ 4 we find that the average number 
of strongly superconducting neighbors of a strongly su¬ 
perconducting site is in the range 0.7 ~ 0.8. Here a 
strongly superconducting site is defined as a site with a 
local parameter Ai > OAAmax where Amax is the largest 
value of Ai in the system (which is close to the maximal 
value Amax = |C^|/2) . Examples of dimers can be seen 
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in Fig. at sites (1,6) — (1,7), (12,15) — (13,15) and 
(8,15) -^9, 15). One also observes that dimers can act as 
seeds of more exte nded i slands as in sites (14, 3) — (14,4). 

In previous workll^l^it has been found that the prefer¬ 
able sites for the SC islands are those with the Hartree 
potential Hi = —\U\{ni) + Vi being close to the chemi¬ 
cal potential /i, since this allows for strong particle-hole 
mixing. This would imply that the ’good’ SC sites are 
already encoded in the normal state since there exists a 
strong correlation between the local Hartree potentials 
in the SC and normal state. On the other hand the cor¬ 
relation between Hi and the size of weakens with in¬ 
creasing disorder, i.e. a small \Hi — /r| not necessarily 
correlates with a large whereas a large A^ always im¬ 
plies a small \Hi — A similar conclusion has been 
drawn in Ref. [26] where the relation of order parameter 
variations and the shell effect has been investigated. 


2. Real space structure of responses 

The largest contribution to the density and amplitude 
correlations comes from the diagonal elements 
Xii that are shown as a logarithmic map in Fig. ^ and 
b, respectively. Here the disorder realization is the same 
as in Fig. and the local SC gap is shown with circles, 
whose size is proportional to the gap magnitude. Panel 
(a) shows also the nearest-neighbor density-density cor¬ 
relation Y(’f encoded in the size of the bars on the bonds. 

One finds that the strong superconducting sites coin¬ 
cide with sites which have a large charge density sus¬ 
ceptibility. Also the dominant nearest neighbor density 
correlations are attached to the SC islands and be¬ 
come particularly enhanced among the sites forming a 
SC dimer. We find that the bare local density correla¬ 
tions Xu^^ ill SC state show a similar structure (not 
shown) but with smaller absolute value (~ 1/20). 

This rises the “chicken and egg” question if sites are 
favorable for superconductivity because they have a large 
susceptibility already in the normal state or if the large 
susceptibility is due to the local superconducting corre¬ 
lations. To answer this question we have computed the 
charge density susceptibility in the absence of supercon¬ 
ductivity. Although there is a tendency for sites with 
charge density susceptibility larger than the average in 
the normal state to become superconducting, there is an 
enormous enhancement of the charge density susceptibil¬ 
ity on the superconducting sites. This can be seen in 
the cut of the local susceptibilities and order parameter 
shown in Fig. We see that on the superconducting 
sites the local susceptibility can be enhanced by two or¬ 
ders of magnitude. The inset shows a zoom of the inten¬ 
sity scale showing that the superconducting sites tend to 
have a charge density susceptibility larger than the aver¬ 
age in the normal state but which does not explain the 
enhancement seen in the superconducting state. It also 
shows that on the sites with small order parameter the 
charge susceptibility remains the same in the supercon- 
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FIG. 3: (Color online) (a): Distribution of the superconduct¬ 
ing gap parameter Ai (circles), the local density correlation 
function Xit (squares), and nearest-neighbor density correla¬ 
tions X^ij) (bars on the bonds), (b): The distribution of the 
local (squares) and nearest-neighbor xftf) (bars on the 
bonds) amplitude correlations together with the SC gap (cir¬ 
cles). (c): Magnitude of local off-diagonal amplitude-density 
correlations \ Xif^\ (squares) together with the SC gap (circles). 
The disorder configuration and parameters are the same as in 
Fig. HI The symbol size for the correlations is displayed on 
a logarithmic scale whereas the SC gap is plotted on a linear 
scale. 
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FIG. 4: (Color online) Cut of the order parameter distribu¬ 
tion (full line, circles) and local density susceptibility x'ii in 
the superconducting state (dashed line, squares) and in the 
normal state (dot dashed line, diamonds). The cut is done 
along the row with t/ = 15 of Fig. and is indicated by an 
arrow in Fig. 


ducting and normal state. Clearly this behavior is due 
to the almost incompressible character of the phase with¬ 
out superconducting correlations which becomes instead 
highly compressible in the superconducting state. This 
physics is similar to that in the clean half-filled Hub¬ 
bard model where a rotation between the two compet¬ 
ing states, CDW and SC, essentially induces a transition 
from zero to very large compressibility k. 

The correlation between SC gap and local charge den¬ 
sity susceptibility is summarized in Fig. which shows 
the distribution of (Ai,Xii) points from 200 samples for 
the normal and SC state and two values of disorder at 
\U\/t = 2. Here always refers to the value in the SC 
state whereas Xu is evaluated in both normal and SC 
state. In the normal state and for weak disorder V/t = \ 
one observes a positive correlation between the local Xu 
and the gap which would develop in the SC state. 
This correlation gets sharper in the SC state (panel b) 
but extends over the same range of Xu values than in 
the normal state. In contrast, for larger disorder Vjt = i 
there is almost no correlation between local charge den¬ 
sity susceptibility and SC gap in the normal state while 
this correlation is strongly enhanced in the SC state and 
pushed to values of yff which are one order of magnitude 
larger than in the normal state. 

The behavior of the amplitude fluctuations is also very 
interesting. We find that local amplitude fluctuations are 
significantly enhanced when the SC gap displays strong 
variations as a function of disorder strength. This feature 
is exemplified in Fig. which, for fixed disorder realiza¬ 
tion (the same as used in Fig. [^and Fig. |^, shows the 
dependence of xti^ on Vb for selected sites. One basically 
observes two kinds of behavior. First there are ’weak’ SC 
sites, as (1,1) or (10, 5), whose order parameter immedi¬ 
ately decreases with the onset of disorder. Besides there 
are ’strong’ SC sites, as (3,12) or (12,15) which initially 





FIG. 5: (Color online) Plot of the points (Ai, xli) (green) for 
the normal (left panels a, c) and SC (right panels b, d) system 
where Ai refers to the value in the SC state. The lines and 
errorbars have been obtained by collecting data in 10 bins of 
A. \U\/t = 2, V/t = 1 (upper panels a, b), V/t = 3 (lower 
panels c, d). 


resist disorder and where A^ can even get enhanced with 
respect to its Vb = 0 value. The drop of A^ on the strong 
SC sites at a given Vb/t is then accompanied by a peak 
in xti^ resembling the behavior close to a second order 
phase transition. However, the order parameter does not 
vanish on the disordered site of the transition but ac¬ 
quires a small finite value due to the proximity effect of 
other SC islands. For a given disorder strength only few 
sites are close to this regime and their number decreases 
with increasing Vo/t due to the decrease of SC islands. 
There are also few sites, as (10,5), where the SC or¬ 
der parameter reemerges at a large value of the disorder 
strength and stays finite over some range of Vb. In the 
appendix the behavior of A^ and xti^ for all sites of the 
sample is analyzed in more detail. 

The present real space analysis reveals that in the 
strongly disordered regime, density correlations are dom¬ 
inant on the SC islands whereas the amplitude correla¬ 
tions are large in the other part of the system, i.e. where 
the SC gap is almost completely suppressed by disorder. 
As shown in Fig. there are only few sites with sig- 
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FIG. 6 : (Color online) Disorder dependence of the SC gap 
(solid, black) and of the local amplitude correlations (red, 
dashed) for selected sites of the disorder configuration used 
in Figs.[^ |f/|/f = 5. 


nificant off-diagonal correlations Xii^ ■ Besides on the 
’marginal’ sites (3,12) and (6,5), which are at the tran¬ 
sition A —>■ 0, the mixing of amplitude and density corre¬ 
lations is only observed on some of the SC sites. Clearly 
this decoupling of amplitude and density correlations will 
be even more pronounced in the average momentum de¬ 
pendent correlations which will be analyzed in the next 
subsections. 


C. Disordered system: Fourier space analysis 


For a particular disorder configuration the Fourier 
transform of the correlation functions is given by 






(17) 


where N denotes the number of lattice sites. Clearly, if 
Xij only depends on the distance between lattice sites 
Ri — Rj then x(q)q0 is diagonal in momenta. In the 
following we perform averages of Xij over different dis¬ 
order realizations up to Ud = 200 for lattice sizes up 
to 24 X 24. This procedure restores translational invari¬ 
ance in the correlation functions so that (x(q, qO)con/. = 
i5(q, q')x(q). In Figs. [8 10 the errorbars in the compress¬ 
ibility and mass reflect the variance of x(q) at q = 0 and 
q = Q, respectively. Although it increases with disorder 
and \U\/t the mean-values exceed the variances for the 
’worst’ cases by a factor ^ 3. 

We fit the correlation function x(q)i which is peaked 
at q = Q, to the function 


^ 1 + 2Ai7i(q — Q)-|-2A272(q — Q) 

with 


(18) 


7 i(q) = 2 - cos{qx) - cos{qy) 

72(q) = 1 - COs{qa:) COs{qy) . 


Although Eq. yields a good account of the cor¬ 

relations over the whole Brillouine zone (BZ) the fit is 
restricted to an area of « 5% of the BZ around the peak 
at Q in order to extract the parameters in Eqs. ( [T5]T^ . 
Expanding Eq. (18) around Q yields 


m 


2 




1 

Aq + A3 

Ai -l- A2 
(Aq -b A 3 ) 
c/ m? = Ai -|- A 2 . 


(19) 

( 20 ) 
( 21 ) 


In the following we analyze the momentum structure 
of the averaged density-, off-diagonal and amplitude cor¬ 
relations. The various fitting parameters will be distin¬ 
guished by (a) the reference momentum in the expansion, 
i.e. q = 0 or Q = (7r,7r), and (b) a superscript which in¬ 
dicates the correlation function. For example, will 
denote the correlation length for amplitude fluctuations 
derived from an expansion of x"^^(q) around q = 0. 



FIG. 7: (Color online) Average (number of samples = 200) 
of Fourier transformed density correlations for parameter 
\U\/t = 2, Vo/t = 3. 


1. Momentum structure o/x^^(q) 

We start with the analysis of the momentum depen¬ 
dence of the averaged density correlation function which 
is shown in Fig.j^for parameters \U\/t = 2 and Vo/t = 3. 
Disorder induces an overall suppression of the response 
as compared to the clean case in Fig. [^. This is most 
pronounced for the CDW correlations at q = Q which 
for Vo/t = 3 are reduced by a factor 1/20 with respect 
to the clean case correlations. At q = 0 this reduction 
is only 1/2 so that in Fig. [^one observes a relative en¬ 
hancement of the zone center correlations. For \U\/t = 2 
the crossover from dominant CDW to q = 0 correlations 
occurs at E/t « 4 whereas for larger values {\U\/t = 5) 
X^'’(q) has a minimum at q = 0 up to the largest dis¬ 
order investigated. Note that also for smaller filling dis- 




















order shifts the dominant correlations from incommen¬ 
surate momenta in the clean case to Q = (tt, tt) so that 
the following analysis is representative for a wide doping 
range and disorder values. 

Fig. ^ shows the parameters (toq)^ and Cq obtained 
from the fit to Eq. 18 with Q = (tt, tt) as a function of 
disorder together with the compressibility k = Xq^o- 



FIG. 8: (Color online) Disorder dependence of the fit parame¬ 
ters (circles), Cq (squares), and (diamonds) for the 

staggered density correlations extracted from Eqs. (18 - 211. 
The disorder dependence of the compressibility is shown by 
the triangles. The dashed-dotted line in panel d) indicates 
the correlation length in the normal state. In panel b) the 
normal state ^g is numerically identical to the result in the 
SC state. 


In the strong coupling limit (small 2t/|17nthe clean 
case compressibility scales as k « \U\/dit^W^ The en¬ 
hancement of K with \U\/t can also be observed in Fig. 
for Vo/t = 0 although the parameters \U\/t = 2,5 are 
rather in the intermediate coupling regime so that the 
agreement with the above estimate is only qualitative. 
Upon increasing Vo/t there is first a decrease of k, in 
agreement with the results of Refs. II9I2I1 At large dis¬ 
order one observes a tendency of the average compress¬ 
ibility K to saturate to a value that is weakly dependent 
on U. Since in this regime the dominant contribution to 
K comes from the (real space) diagonal elements Xii on 
the SC islands there exists an apparent inverse correla¬ 
tion between the number of SC islands (which decreases 
with Vo/t) and the local compressibility Xii (which gets 
enhanced with increasing Vo/t). 

We now turn to the analysis of the CDW correlation 
length in the disordered SC system. For weak disorder 
Vo/t = 0.5 there is a strong difference in the density 
distribution obtained for the two values of \U\/t = 2,5 
which we have investigated. In fact, for \U\/t = 2 we 
find that the difference in the density distribution be¬ 
tween normal and SC state is small for each value of the 
disorder potential Vo/t. As a consequence the decrease 
of the CDW correlation length with Vo/t (Fig. [^) is the 


same in the normal and SC state within the numerical 
accuracy. On the other hand, for \U\/t = 5 we find that 
already for Vo/t = 0.5 sites in the normal state system 
are either almost empty or doubly occupied. As already 
discussed above, the SC state induces a redistribution of 
charge density which in this case leads to a significant re¬ 
arrangement with a more homogeneous distribution be¬ 
tween n « 0.2 and n « 1.7. As a consequence of this 
effectively less disordered SC state one observes in panel 
(d) of Fig. [^an enhancement of the correlation length at 
Vo/t = 0.5 from « 0.3 in the normal state to ^g « I 
in the SC system. 

The behavior of fit parameters in the SC system, 
as shown in Fig. can then be qualitatively under¬ 
stood from the evolution toward the bimodal charge den¬ 
sity distribution, where the low (high) density peak ap¬ 
proaches riL = 0 (uh = 2) with increasing disorder. We 
also adopt the result from a strong coupli ng expansion 
of foi' homogeneous systerr(27li2l which for the 

mass parameter yields 


2 8t2 ,52 


( 22 ) 


and 5 = 1 — n denotes the doping measured from half¬ 
filling. Averaging Eq. over the bimodal distribution, 
yields (wg) ~ {5n)‘^/{\ — (^n)^) with Sn = — ul- 

The grow of 5n with Vo/t then accounts for the increase 
of TOg with disorder as shown in Fig. 

In the strong-coupling clean case the parameter cq is 
given bjE^ 


_ 1 _ 2(5^ ^2 

= |[ 7 | 1-52 = |; 7 | “ ^ 


(23) 


and is thus expected to decrease with disorder propor¬ 
tional to the increase of uiq . Within the numerical error 
this is in fact the behavior observed in Fig.|^and also ac¬ 
counts for the decrease of the correlation length with 
disorder. 


2. Momentum structure of amplitude correlations 

We proceed by analyzing the amplitude correlations 
^AA 

(q) on top of the BdG solution whose momentum 
dependence is reported in Fig. |^for|C/|A = 2, UoA = 3. 

It turns out that disorder removes the enhancement of 
amplitude correlations at Q = (7r,7r), which were domi¬ 
nating in the clean case for this value of \U\/t. An inter¬ 
esting result is the concomitant enhancement of the q = 0 
response by a factor of ~ 5/2 which therefore dominates 
the amplitude correlations for large disorder. As we have 
seen in the previous section, the density correlations are 
still peaked at Q = (tt, tt) for these parameters which 
indicates the decoupling of density and amplitude fluc¬ 
tuations with increasing disorder. Note that in contrast 
to the density correlations, the amplitude fluctuations in 
the normal state will always be unstable. 
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FIG. 9: (Color online) Average (number of samples = 200) 
of Fourier transformed amplitude correlations for parameter 
\U\/t = 2, Vo/t = 3. 


The latter are again characterized by the mass (m^) 
and cjf parameter obtained from the fit of fo 

Eq. 18 around q = (0, 0). Fig. [^reports the fit param¬ 
eters as a function of disorder, again for values of the 
onsite attraction \U\/t = 2 and \U\/t = 5. Note that 
for the larger interaction \U\/t = 5 and small disorder 
the correlations show the dominant peak at Q = (tt, tt) 
for which reason the fit parameters are only reported for 
Vo/t > 0.5. 


The aforementioned enhancement of the q = (0,0) am¬ 
plitude correlations with Vb /t now results in the decrease 
of the mass with disorder with tendency to satu¬ 
rate at large Vo/t > 2. Also the parameter c decreases 
with the disorder strength so that the resulting correla¬ 
tion length = Co/{rnoY (right insets to Fig. [M cru¬ 
cially depends on the relative change of Cq and(m;^)^ 
with Vo/t. 

For \U\/t = 2 the correlation length is almost constant 
up to Vo/t = 2.5 and then starts to decrease with disor¬ 
der. For larger \U\/t one even observes an enhancement 
for small Vo/t so that ^o acquires a maximum around 
Vo/t = 2.5. We note that this is not an effect of compet¬ 
ing CDW order since the same result is observed in the 
low-density regime where such correlations are absent. 

In the limit of small Vo/t one can adopt the usual 
expression for the correlation length in dirty supercon¬ 
ductors given by = VCbcsI with the mean free 
path I and the correlation length of the clean system 
^BCS ~ vp/^^^■ The behavior of ^o(^) therefore cru¬ 
cially depends on the depletion of the density of states, 
which lowers the superconducting gap, and the re¬ 
duction of the mean free path I with disorder. As noted in 
Ref. ll9l2T] the situation in the strongly disordered system 
is more interesting since one has to distinguish between 
the average superconducting order parameter and 

the spectral gap. As shown in the left insets to Fig. [I^ 
(A^C) continuously decreases with disorder due to the 
increase of the ’non-SC’ area. On the other hand the 
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FIG. 10: (Golor online) Disorder dependence of the fit param¬ 
eters (mO^)^ (circles), Cq (sqnares) and = Y(Y/(m-^Y 
(diamonds, right inset) as extracted from Eqs. (18-211 for 


\U\/t = 2 (a) and \ U\/t = 5 (b). The left inset reports the av¬ 
erage superconducing gap (circles) and average spectral gap 
(squares). The right insets also show the gap autocorrelation 
length Aac (circles) computed from Eqs. ( |24|^ . 


spectral gap first shrinks with disorder due to the deple¬ 
tion of the density of states but grows again for strong 
disorder, signaling the formation of local boson pairs that 
get progressively localised as the SIT is approached. One 
can then argue that at strong disorder the BCS correla¬ 
tion length tends to scale as the inverse of the spectral 
gap, that acts as a cut-off to the increase of associated 
to the suppression of the SC order parameter. Alterna¬ 
tively one can relate the disorder dependence of the cor¬ 
relation length to the behavior of the nearest-neighbor 
amplitude correlations as shown in the appendix. 

Recentljff^ the spatial dependence of the STM spectra 
in strongly disordered NbN films has been analysed in 
terms of the autocorrelation function for the order pa¬ 
rameter, i.e. 

- (^)) - (^)))- ( 24 ) 

By performing an average over several disorder configura- 
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FIG. 11: (Color online) Average of Fourier transformed off- 
diagonal correlations for = 2.0 and |17|/t = 2. 


tions we can extract the corresponding correlation length 
Xac from a ht to the function 

F(R) = ao + sin^(2<i>)-ta3 sin=( 40 )) ^ (25) 

Here (j) is the polar angle related to R which incorpo¬ 
rates anisotropies in the correlations and we restrict the 
fit to |R| > 2 in order to isolate the long-distance behav¬ 
ior. The resulting length Xac as a function of disorder is 
shown by circles in the right inset to Fig.j^and it is close 
to the correlation length extracted from the amplitude 
correlations. For Vq —t 0 one can apply linear response 
theory on the disorder and show that the two lengths co¬ 
incide. In the strongly disordered regime the situation is 
more complex. We hnd numerically that both lengths are 
close to each other. Notice that for \U\/t = 5 we observe 
that both Xac and increase in the regime where the 
separation between the order parameter and the spectral 
gap starts to develop, while they collapse in the regime 
where the spectral gap tends to increase again. In Monte 
Carlo simulationsSl the latter regime corresponds to the 
SIT, not captured by the present Bogoliubov-de-Gennes 
approach. This same tendency is observed in the experi¬ 
mental estimate of Xac given in Ref. [T^l done for samples 
in the so-called ’’pseudogap” region of the phase diagram, 
where the spectral gap is much larger than Tc- 


3. Momentum structure of off-diagonal correlations 

Off-diagonal correlations mix the density and 

amplitude sector and are shown in Fig. for the clean 
case and in Fig. [^for the disordered system. 

Upon coupling an external field in the density sector 
Hi — J2ci\P-q. correlation function yields 

the corresponding response for the gap amplitude. In 
particular, for q = 0 a spatially constant (and positive) 
Aq=o induces an effective reduction of the chemical po¬ 
tential. Consider now the clean case where for the attrac¬ 
tive Hubbard model with nearest-neighbor hopping the 


gap amplitude as a function of density has a maximum at 
half-filling and continuously decreases towards n = 0 and 
n = 2. Therefore off-diagonal correlations are negative 
for n < 1 (where a positive A shifts the effective chem¬ 
ical potential away from half-filling) in agreement with 
Fig.[T] and positive for n > 1. Similar arguments can be 
made for finite momenta. In particular, the strong en¬ 
hancement of |x"^^(q)l at q = Qcdw observed in Fig. 
is due to the strong competition between CDW and SC 
correlations close to half-filling. 

In the doped system Fig. E] reveals a strong sup¬ 
pression for the off-diagonal correlations due to the spa¬ 
tial separation of density- and amplitude fluctuations as 
demonstrated in Sec. |HI B| Naturally this is again most 
pronounced for the CDW momentum due to the removal 
of particle-hole symmetry by disorder. It is worth noting 
that in the dynamic limit (q = 0, a; finite) the off-diagonal 
correlations show instead the opposite behavior. More 
specifically, as it has been recently discussed in Ref. [33] , 
the coupling between the amplitude and density/phase 
correlation at finite frequency is strongly enhanced by 
disorder, leading to a strong mixing between the ampli¬ 
tude and phase spectral functions at zero momentum. 


IV. CURRENT CORRELATIONS 


To conclude our analysis of the SC correlations we shall 
discuss now the change in the current-current correla¬ 
tion function induced upon entering the superconducting 
state. In particular we want to explore the consequences 
of the percolative current formation (cf. Fig. on the 
behaviour of the current correlation function entering 
the definition (12) of the superfluid stiffness. 


In order to obtain the intrinsic superconducting re¬ 
sponse, we have to subtract the contribution which is 
already present in the normal state (at finite momenta), 
and which can be either diamagnetic or paramagnetic 
depending on the filling of the system. 

This is illustrated by the dashed line marked 
with diamonds in Fig. for the homogeneous non¬ 
superconducting system. Clearly, the current response 
of Eq. ([^, Dq^ = -{Tafj -\- {x^^Iqv)) vanishes at = 0 
(i.e. the SC stiffness) when the system is in the normal 
state, however, it becomes non-zero for hnite momenta. 


In particular at low density (cf. Fig. 12 1 ) one recovers 
the finite-q diamagnetic response {Dq^ > 0) related to 
Landau diamagnetism in agreement with the transverse 
current response of a Fermi liquid.^ In contrast, larger 
filling (cf. Fig. |T^) supports a finite-q paramagnetic 
current response which would even diverge at q = (tt, tt) 
for n = 1 (not shown). This feature is the starting point 
for the exploration of circulating current phases as pos¬ 
sible candidates for the pseudogap in cuprate supercon- 
ductors.l^ 

As shown by the triangle symbols in Fig.[T^a finite SC 
gap shifts up the curves in order to yield a diamagnetic 
Dq^ independently on doping. In order to extract what is 
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FIG. 12: (Color online) Transverse current response Dq^ = 
~{Tx} + {x^^ {%)} for the non-SC (blue dashed, diamonds) and 
the sc homogeneous system at \U\/t = 5 (red dashed, trian¬ 
gles). The normal state response (Asc = 0) is independent of 
\U\/t. The solid lines report the difference ADa{qy) between 
Dq^ for the sc- and normal system for \U\/t = 5 (squares) 
and \U\/t = 2 (circles). Filling n = 0.325 (a) and n = 0.875 
(b). 


due to superconductivity we take the difference with re¬ 
spect to the normal state response corre¬ 

sponding curves are shown by square symbols (|t/|/t = 5) 
and circles {\U\/t = 2) in Fig. 12 for n = 0.325 and 
n = 0.875, respectively. In the weak coupling limit the 
difference ADs{qy) = jg ^^.^yayg strongly 

peaked at = 0 and the underlying normal state re¬ 
sponse does not influence the curvature of the peak which 
determines the SC coherence length. On the other hand, 
it turns out that for large filling and strong coupling (cf. 
squares in Fig. 121 ADs{qy) can even acquire a maximum 
at the zone boundary. Thus in this limit the SC diamag¬ 
netic response is largest on short length scales and corre¬ 
sponds to an oscillatory decay of the SC induced current 
correlations in real space. 

shows the transverse current response ADs{qy) 


Fig. 


13 


for various disorder strength and interaction \U\/t = 2 
with the normal state result substracted. The latter has 
obtained for the same disorder configurations and by set¬ 
ting Af^ = 0. 



‘Jy 


FIG. 13: (Color online) Main panel: Transverse current corre¬ 
lations ADs{qy) measured with respect to the normal system 
for \U\/t = 2, n = 0.875, and various disorder strengths. Up¬ 
per left inset: Correlation length extracted from Eq. (26l. 
Upper right inset: Superfluid stiffness. 


We parametrize the long-wavelength structure as 

ADa{qy)=Da[l-iCDqyf] (26) 

which defines a SC coherence length related to the dia¬ 
magnetic response and allows us to extract the stiffness 
Da as a function of disorder. Both quantities are shown 
in the insets to Fig. 

As discussed previousljEH (see also Sec. Ill B 1) Ds gets 


rapidly suppressed with disorder but since the BdG ap¬ 
proach does not capture the SC-insulator transition it 
does not vanish even for large Vo/t. Also the coherence 
length (cf. left inset to Fig. 13) is strongly suppressed 
by disorder. Above Vo/t « 2/AD{qy) is essentially in¬ 
dependent on the transverse momentum qy and « 0 
within the numerical accuracy. 

However, due to the average over disorder configura¬ 
tions the above analysis does not capture the long-range 
current correlations which exist along the percolative 
path (cf. Sec. IIIBI and which we will analyze separately 
in the following. 

First we identify the superconducting backbone. The 
criterion to decide which sites belong to the percolative 
path is chosen as follows: For the vector potential A 
along the x direction, we determine the maximum cur¬ 
rent through a bond in the system and select all 

sites which have currents larger than We find 

that usually a value of a = 1/3 is appropriate in order to 
selecting the sites which are visited by the path. An ex¬ 
ample is shown in the inset to Fig. |14| where the squares 
indicate sites with jx{Rn) > /^- Clearly, there are 

sites (e.g. in the upper right corner) which are traversed 
by a minor current but are left out by the ’a = 1/3’ crite¬ 
rion. Reducing further the value of a would also include 
these sites, however, we note that the following results do 
not depend sensitively on the value of a. The effect of a 
larger (smaller) a is to add sites with larger (smaller) cur¬ 
rent to the path which concomitantly slightly increases 
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(decreases) the long-distance correlations which are cal¬ 
culated below. 
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FIG. 15: (Color online) Current correlations {ADnm) aver¬ 
aged over 200 samples for disorder strength V/t = 3 and 
n = 0.875. Full symbols: |?7|/t = 2; Open symbols: |17|/t = 5. 


FIG. 14: (Color online) Main panel: (ADnm) for both sites 
(black) and only one R„ (blue) on the percolative path shown 
in the inset. The squares indicate sites R„ with jx{Rn) > 
|17|/t = 2, n = 0.875, V/t = 3. 


by using the following classical phase-only action: 

(28) 

i 


We proceed by evaluating the non-local stiffness Dn^m 
between sites i?„ and Rm 

DTm = [-Sn,mtx{n) - Xnm{j/l,j/n)] ( 27 ) 

and compute the difference between sc and normal state 
AHrim = 71®^ — Two cases are considered: (a) 

both sites i?„ and Rm belong to the percolative path 
and (b) only one of the sites i?„, Rm is on the path. 
The result for D^m in both cases is shown in Fig. 
for the particular percolative path displayed in the inset. 
The ’errorbars’ indicate the variance due to the fact that 
different sites i?„ and Rm have the same distance |i?„ — 
Rm\ but different values for D^m- 

As can be seen the current correlations rapidly decay 
away from the percolative path and are practically ’zero’ 
for |i?„ — Rm\ > 3. On the other hand correlations on 
the path stay finite up to the largest distances available 
in the system. 

Finally, Fig. [TS] shows the on- and off-path current cor¬ 
relations averaged over 200 disorder configurations for 
|17|/t = 2 and \U\/t = 5, respectively. As for the spe¬ 
cific sample shown in Fig. the off-path correlations 
get rapidly suppressed while on-path correlations stay fi¬ 
nite up to large |i?„ — Rm\- Upon comparing the on-path 
correlations between the two |17|/t-values one finds, be¬ 
sides a reduction by a factor of « 10, that the decay of 
ADnm with distance for |t7|/t = 5 is significantly smaller 
than for |17|/t = 2 while still staying finite for the largest 
possible separation in the system (« y/2 x 10 for a 20 x 20 
lattice when the percolative path is along the diagonal). 

The persistance of the current correlations along the 
percolative path resembles closely the expected behavior 
for a one-dimensional chain, where it simply follows from 
the current conservation. This can be easily seen at w = 0 


where Ji are the local (random) stiffnesses (in units of 
the temperature T) and represents the local phase 
gradient (5$^ = (6>i+i — 0i), 9i being the local SC phase. 
Eq. (28) can be obtained for example by expanding at 
Gaussian level a classical XY model with random cou¬ 
plings Ji, that is the prototype model for the phase de¬ 
grees of freedom of a superconductor. Eq. (28) is also 
obtainecPS by mappin^2ll at large U the disordered Hub¬ 
bard model into the pseudospin model. In this mapping 
the superconductivity corresponds to a spontaneous in¬ 
plane magnetization, i.e. to the usual XY model with 
a coupling J ~ /U, and disorder maps into a random 
out-of-plane field, that leads in turn to the disorder in 
the local couplings Ji after a Holstein-Primakoff expan¬ 
sion around the mean-field solution.l^ 

The local current li for the model ( |2^ can be written, 
after minimal coupling substitution — 2Ai in 

Eq. (1^ as: 


h = 2J,{5^, -2A,). (29) 


In the one-dimensional case the current conservation 
implies that li is independent on the site index, i.e. 
(i5d)i — 2Ai) = c/2Ji, where c is a constant. By sum¬ 
ming over the site index and using the boundary condi¬ 
tion = Oone then gets c = -'^{Y.i^i)/Ji)- 

Since the superfluid stiffness is defines as usual (see Eq. 
([I4|)) as Ds = —I{q = 0)/A{q = 0) one also deduces that 



so that li = c = —{^/N)J2jDsAj. By comparing this 
with Eq. (27) above we then recover that Dij = Dg/N 
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for all pairs of sites i,j along the chain. It is interest¬ 
ing to note that this result also implies that the param¬ 
agnetic contribution to the current must cancel out the 
local diamagnetic term AJi of Eq. (29). This can be seen 
by computing explicitly the average current value from 
Eq. (29) in linear response theory, in analogy with the 
expression ([27|) introduced above: 


(7j) — 4 Jj{Sjj XijJj)Aj — Djj Aj (31) 

j 3 


where Xij = {5^i5^j) is easily determined from Eq. 
as: 


(28) 


{S<A>,6<A>,) = 

f dAVS^exp [-4 Jk{^^kY + ^^k] 

/ dXVS^ exp [-X Jk{5^kY + *A Y.k 

(32) 


where the A integration accounts for the periodicity con¬ 
straint. By making the change of variables (5$^ —>■ 
— iX/ Jk one immediately sees that 


Y _ ^j3_ _ 

~ J, NJ^Jj 


(33) 


that inserted into Eq. (31) gives Dij = Ds/N, as antici¬ 
pated before. We note mso that the independence of Dij 
in Eq. (311 on both site indexes can be also derived as a 


consequence of charge conservation and gauge invariance 
in one dimension. Indeed the independence of Dij on the 
site index i is a consequence of a constant current li on 
each site, while the independence of Dij on the second 
index j is a consequence of the fact that at w = 0 only 
to the g = 0 component of the gauge field A leads to a 
finite response. 

Going back to our 2D system, we clearly see in Fig. [I5| 
that for a fixed disorder strength the percolative path 
becomes more ’ID’-like with increasing \U\/t, which ac¬ 
counts for the crossover to a more constant AD^m for 
\U\/t = 5. Indeed, a larger \U\/t corresponds to a smaller 
J in the mapping into the XE-like bosonic model, with 
an enhanced influence of disorder and with smaller ef¬ 
fective local stiffnesses J^. This in turn is in agreement 
with the strong reduction of ADnm from \U\/t = 2 to 
\U\/t = 5, as shown in Fig. 15 


is observed in disordered films of conventional supercon¬ 
ductors, like e.g. NbN, InOa, and TiN.lsESllt is then cru¬ 
cial to assess how this inhomogeneous SC state affects 
the behavior of the amplitude, density and current cor¬ 
relations, in order to interpret the results of the various 
experimental probes. 

In the present manuscript we analyzed this issue within 
the fermionic Hubbard model with on-site disorder. We 
presented a detailed study of the correlation functions 
both in real space, for a specific disorder configuration, 
and in momentum space, after the average over several 
disorder configurations. The momentum-space analysis 
allows us to extract the correlation length of each physi¬ 
cal quantity in close analogy with the usual approach for 
homogeneous systems. As a first result, one then sees 
that while in the homogeneous case at low temperature 
amplitude and current correlation lengths coincide up to 
a numerical factoJ^, in the presence of strong disorder 
this is no more true as can be seen in the summariz¬ 
ing figure By means of a simultaneous analysis of 
the real-space correlations we can then disentangle how 
the properties of the fragmented SC ground state influ¬ 
ence the various correlation lengths. As we discussed in 
the manuscript, these two approaches give complemen¬ 
tary informations, that we will summarize below. In this 
respect, even though our results are based on a RPA ap¬ 
proximation, they have the advantage to allow for larger 
system sizes than Monte Carlo simulations, as e.g. those 
reported in Ref. [Ml The use of large clusters is in turn 
crucial to trace back the behavior of different response 
functions to the inhomogeneous structure of the ground 
state and to perform a momentum-space analysis. 



V. DISCUSSION AND CONCLUSIONS 


As we discussed in the introduction, it has been now 
established in several theoretical models that when the 
SIT is approached a granular SC state emerges, with SC 
puddles embedded in a non-SC background. Thanks to 
the enormous progresses made in the experimental tech¬ 
niques able to probe the systems in real space, it has 
been also established that such an emergent granularity 


FIG. 16: (Color online) Summary of results for the various 
correlation lengths as a function of disorder for \U\/t = 2 
and n = 0.875. At very small disorder the autocorrelation 
length Xac cannot be properly defined since the approximated 
formula ( |25| ) does not reproduce accurately the data, see also 
C{R) in Fig. 17 below. 


Amplitude and density correlations. 

We find that in general the strength of the amplitude 
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response ~ l/(m^)^ increases with disorder while the 
charge response ~ l/irriQ)'^ gets suppressed by disorder 
(cf. Figs. [^ari d[l0|). This is similar to a previous Monte 
Carlo study^^ which found that superconducting correla¬ 
tions are much more robust to disorder than charge cor¬ 
relations. Here, due to the larger system size, we could 
explore in detail the origin of this behavior. 

The suppression of the charge response is easily un¬ 
derstood by the tendency of disorder to localize the pairs 
and render the system incompressible almost everywhere 
except in the superconducting islands. The increase 
of the superconducting response is more subtle. For 
strong disorder the region in between the islands contains 
“marginal” sites where the order parameter is small but 
very susceptible to become large by small variations of 
the disorder [see Fig.j^b) and Fig. I^for site (3,12)] yield¬ 
ing a large overall pair susceptibility and resembling the 
behavior close to a second order phase transition. The 
decoupling of density and amplitude correlations in real 
space is reflected in the momentum-space structure of the 
susceptibilities. Thus, while in the homogeneous cas^^ 
the maximum or Xpp(q) at the CDW vector Q = (7r,7r) 
leads to an enhancement of the amplitude correlations 
XAAi^) at the same wavevector (see Fig. [^, in the dis¬ 
ordered case this effect disappears (Fig. [^. 

The resulting amplitude correlation length shown 
in Figs. has an interesting disorder dependence. 

Indeed, it stays constant or it is even enhanced at in¬ 
termediate disorder levels, before then being ultimately 
suppressed as the SIT is approached. In the latter regime 
we argued that the decay of the correlation length is ruled 
by the behavior of the spectral gap, which increases as 
pairs become localized with disorder. 

In Figs. 10 we also compared with the autocorrela¬ 
tion length Xac, that can be directly extracted experimen¬ 
tally from the STM maps of the SC ground state. This 
has recently been done for disordered NbN film^i^ and we 
show for convenience the corresponding data in Fig. [T^l. 
In this work the SC islands are identified by the regions 
with a lar ge SC coherence-peak height, that is usually 
takeiP^^^ to be a measure of the local order parameter 
Ai, i.e. the local gap solution in the BdG equations. By 
analyzing the spatial correlations between good SC sites 
the authors of Ref. [T2] found that the autocorrelation 
length Xac becomes larger as disorder is increased. This 
is shown in Fig. where we report the experimental 
data for the autocorrelation function C{R) defined in Eq. 
(24) above. A similar trend can be observed also in our 
simulations, see Fig. 17 3, where C{R) shows first a rapid 
suppression over a length scale of the order of the SC 
island, followed by a long-tail decay that can be eventu¬ 
ally fitted with the approximated formula (25) in order 
to extract Xac- Since this tail can be thought as the re¬ 
sponse of the system to the fluctuations that created the 
island we expect that Xac is close to as indeed we find 
numerically, see Fig. and Fig. 

In contrast to the autocorrelation length, a direct esti¬ 
mate of the amplitude correlation length from the ex- 



FIG. 17: Comparison between the experimental estimate (left 
panel) of the antocorrelation function, defined in Eq. ( |24[ |, 
and the numerical computations (right panel). The experi¬ 
ments data are taken from Ref. m and refer to three NbN 
films at different disorder level (labeled by the different crit¬ 
ical temperatures Tc). The theoretical data are obtained for 
\U\/t = 2, n = 0.875, and disorder values Vo/t = 0.5 (solid) 
and Vb/t = 3. (dashed). Considering that the typical size of 
the SC islands in these NbN films range between 20-40 nm, 
and it is one-two lattice spacings in our simulations, the length 
scales in the experiments and simulations are approximately 
comparable. 


periments is not so straightforward. Indeed, while within 
a Ginzburg-Landau approach, where a single length scale 
exists, ^0 can be estimated from the upper critical field at 
T = 0 as Hc 2 = ‘ho/(27r^o)> ^-t strong disorder this con¬ 
nection is not obvious. In particular when the superfluid 
stiffness Dg is the lowest energy scale in the problem one 
would expect that Tc oc Dg^ so that also the upper criti¬ 
cal field will scale with Dg, as suggested for example by 
a recent analysis of the microwave conductivity at finite 
magnetic field in disordered InOj, In this sense, even 
though at intermediate disorder the decrease of Hc 2 mea¬ 
sured experimentalljEl can be interpreted as an increase 
of ^0 due to the weakening of the SC order parameter, 
as the SIT is approached one should not attribute the 
vanishing of Hc 2 oc Tc to a divergence of discussed 
above. 


Current correlations The behavior of the current cor¬ 
relations is also strongly influenced by the formation of a 
fragmented SC state. Indeed, as already noticed before,!^ 
the superfluid response is mainly determined by a few 
percolative paths that connect the good SC regions. As 
a consequence, the decay of the current correlations de¬ 
pends on the position of the initial and final sites with re¬ 
spect to this SC ’backbone’. If both sites belong to a per¬ 
colative path the current correlations are long-ranged (es¬ 
sentially constant, see Fig. 14), in agreement with what 
one expects for a truly one-dimensional system, like e.g. 
the one-dimensional XY -model. On the other hand, this 
long-range behavior is easily missed when the transverse 
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current correlations are extracted from the response in 
momentum space after average over several disorder con¬ 
figuration. Indeed, the current-current correlation length 
is rapidly suppressed (cf. inset to Fig. 13 and Fig. 
16 1 ), in analogy with the overall superfluid response. 


This behavior has to be contrasted to the one of the am¬ 
plitude correlation length that is strongly suppressed 
only at the SIT. On the other hand, the persistence of 
current correlation along the percolative paths suggests 
that the existence of the SC backbone can be deduced 
in principle by the measurements of the space-dependent 
current susceptibilities, without having to evaluate ex¬ 
plicitly the current pattern at finite applied field. The ex¬ 
perimental study of these issues is of course challenging, 
but it should be accomplishable with four-point atomic 
force microscopy when the electrode spacing reaches the 
nanometer separation. Its observation would certainly 
contribute signihcantly to our understanding of the basic 
mechanisms leading to the formation of the inhomoge¬ 
neous SC state as the SIT is approached in real systems. 
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Appendix A: Disorder dependence of local SC gap 
and local correlations 


Fig. [^reports the disorder dependence of the SC gap 
value and local amplitude correlations on each site for the 
same disorder configuration and parameters used in Figs. 
Hi Note that the amplitude correlations are normalized 
to their maximum value at each site. Clearly, the SC 
order parameter on the majority of sites drops to a small 
value around Vb/t « 2 but there are also singular sites 
where extends up to Vq/^ ~ 4 or where reemerges 
at large disorder values. 

In the clean system the onset of a hnite SC gap be¬ 
low Tc is accompanied by a divergence in the amplitude 
correlations, both the local and non-local ones. The 
pronounced enhancement of the amplitude correlations 
around Vq/^ ~ 2 in Fig. 18 3 suggests a similar fea¬ 


ture as a function of disorder with the difference that A^ 
does not vanish but becomes small beyond some value 
of Vq- To analyze this feature in more detail we plot 
in Fig. 19 the probability density P(A < e) as a func¬ 
tion of the disorder strength Vq- Here P(A < e)dVo is 
the probability that the order parameter of a given site 
will fall below the threshold e for the first time when 
the disorder is increased from Vq to Vq + dVo. Also 
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FIG. 18: (Color online) Top panel: disorder dependence of 
the SC gap value for the same disorder configuration used in 
Figs. The site index for the 16 x 16 lattice is obtained 

from ix + 16{iy — 1). Lower panel: disorder dependence of the 
local amplitude correlations xfi'^ normalized to their maxi¬ 
mum value at each site. \U\/t = 5, n = 0.875. 


shown are the probability distributions for the maxi¬ 
mum in the local [P{xu^ = max)] and nearest neighbor 
= max)] amplitude correlations where, for ex¬ 
ample, Pixti^ = max)dVo is the probability that Xu^ 
for a given site i, attains its maximum value as a func¬ 
tion of disorder in the interval Vq, Vq + dVo- Clearly, 
for ]U]/t = 5 (right panel of Fig. P(A < O.Olt) 
has a pronounced peak around Vo/t ~ 2 ... 2.5 and one 
finds that for about 50% of all sites A^ < O.Olt between 
1.5 < Vo/t < 2.5. Concomitantly also the probability dis¬ 
tributions for the local and non-local amplitude correla¬ 
tions are peaked at a somewhat lower value of Vq/^ ~ 1-5. 
For smaller ]U]/t = 2 these distributions are broader and 
in particular the nearest-neighbor amplitude correlations 
are no longer characterized by a significant enhancement. 
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This finding offers an alternative perspective for under¬ 
standing the disorder dependence of the amplitude corre¬ 
lation length ^0 shown in Fig. Since is of the order 
of one lattice spacing the nearest-neighbor correlations 
yield the dominant contribution to the correlation length 
which accounts for the enhancement around Vo/t = 2. 
On the other hand, the distributions as a function of 
Vb/t are significantly broader for \U\/t = 2 (cf. Fig. 19 1 ) 
which agrees with the behavior of shown in Fig. 10 r. 


FIG. 19: (Color online) Black (full) steps: Probability distri¬ 
bution P{A < e) that the order parameter of a given site will 
fall below the threshold e for the first time, upon increasing 
the disorder with e = O.Olt ; Red dashed step: Probability 
distribution P{xti^ = max) that the local amplitnde correla¬ 
tion of a given site will attain its maximum value as a function 
of disorder strength. Blue thin step: Probability distribution 
P{xfif) = max) that the nearest-neighbor amplitude corre¬ 
lations of a given bond will attain its maximum value as a 
function of disorder strength. Left panel: \U\/t = 2, Right 
panel: \U\/t = 5. 
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